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We analytically solve the problem of Bose star growth in the bath of gravitationally interacting 
particles. We find that after nucleation of this object the bath is described by a self-similar solution 
of kinetic equation. Together with the conservation laws, this fixes mass evolution of the Bose star. 
Our theory explains, in particular, the slowdown of the star growth at a certain “core-halo” mass, 
but also predicts formation of heavier and lighter objects in magistral dark matter models. The 
developed “adiabatic” approach to self-similarity may be of interest for kinetic theory in general. 


1. Introduction. Gravitationally bound blobs of 
Bose-Einstein condensate [1] — Bose stars — have re- 
gained a lot of attention recently. This is because they 
are abundant in models with light dark matter [2] consist- 
ing, e.g., of “fuzzy” bosons or QCD axions. In those two 
cases, the Bose stars are called “solitonic galaxy cores” [3] 
and “axion stars” [2], respectively. Typically, the self- 
couplings of light dark matter particles are tiny and can 
be ignored. But their phase-space density is so large [4] 
that thermalization can occur inside the smallest cos- 
mological structures via universal gravitational interac- 
tions [5]. This makes the Bose star appear in the center 
of every such structure [3, 5, 6] in kinetic time. 

The question is, how do the newborn Bose stars grow? 
Lattice simulations show that their masses increase at 
first as My, x t'/? [5] and then slow down [6]. But the 
numerical results on the late-time behavior are conflict- 
ing: My, œ t!/® in [6, 7] and ¢1/¢ in [8]. 

In this Letter, we for the first time show [9] that Bose- 
Einstein condensation of dark matter via gravitational 
(long-range) scattering is described by self-similar solu- 
tions of kinetic equation. Computing the condensation 
flux onto the Bose star, we analytically obtain its growth 
law. The star mass is not a simple power of time, but can 
be piecewise approximated by all of the behaviors above. 

2. A crucial observation. Consider a cloud of nonrel- 
ativistic dark matter bosons inside the smallest cosmo- 
logical structure: a minicluster of axions [10] or a galaxy 
halo of “fuzzy” dark matter. At small particle masses m 
the occupation numbers are so large that the bosons are 
described by a random classical field y(t, æ) evolving in 
its own gravitational potential U(t, x), 


ibp = -AY /2m + MUY, (1) 
AU = 47G (mlwl? — P), 


where p = M/L? is the mean density and M is the total 
mass. For simplicity, we approximate the structure with 
periodic box of size L. 

We solve Eqs. (1) using a stable 3D code of Ref. [5]. 
The starting point of this evolution is a virial equilib- 
rium, i.e. Gaussian-distributed field with Fourier im- 
age |p|? a Me~?’/?6 and random phases arg Wp; here 
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FIG. 1. Simulation of Eq. (1) with M = 50 po/m?G and L = 
60/po. (a) Spectra (2) at two moments of time (solid lines). 
(b) Bath spectra (w > 0) at large times and (c) their self- 
similar transformation (3) with D = 2.8. Figure (b) includes 
t = tgr graph (dash-dotted). Chain points in Figs. (a), (c) 
show the solution of Eq. (5) with D = 2.8 and the source; see 
Supplemental Material C (SM-C) for parameters. 


wo = pê /2m is the typical particle energy [11]. We study 
mass distribution F(t, w) = dM/dw of bosons over ener- 
gies w. It is given by time Fourier transform [5]: 


dt’ . $ f 
. =m f z- Pey (t eplet, g) A, (2) 


where At™! < wo is the energy resolution. In the 
isotropic homogeneous case, this function is related to 
the usual phase-space density f(p) as F = L’m?pf/2r° 
with w = p?/2m. Below we exploit dimensionless units: 


F(t, ©) = 2woF/M and ð = w/2wo. 

Our ensemble has large occupation numbers and hence 
thermalizes into a Bose-Einstein condensate. This is seen 
in simulations [5] as phase transition at a kinetic time 
t = tgr = 2bV2 w8 /[3r3 G? In(poL)], where b ~ 0.9 for 
the Gaussian initial distribution. Namely, after tgr the 
spectrum F(t, w) develops a narrow peak at w ~ wy, < 0 
moving with time to lower energies; see Fig. 1(a) and the 
video [12]. The peak is a Bose star [1, 2]: a condensate 
of particles occupying a single — ground — level wy, in 
the collective gravitational well U. Once the Bose star 
appears, the ensemble mass M = Mys + Me + Mep divides 
between this object (Mbs), excited bound states in its 
gravitational field (Me), and the “bath” of particles with 
w > 0 (M,). The conditions for condensation are still 
satisfied, so M, grows at t > tgr. 

Below we measure time in kinetic intervals T = t/tgr 
and compute M; by integrating F'(w) over the respective 
regions. E.g., the “dressed star” mass is M, = My, + 
Me = f <o Fdw, while M, and Mes are the integrals 
over w > 0 and w & wps, respectively. 

Now, we make an important observation. Consider 
the w > 0 spectrum, ie. the “bath”. It changes a lot 
after the Bose star formation, cf. the graphs with differ- 
ent 7 in Fig. 1(b). The same graphs, however, coincide 
in Fig. 1(c) after time-dependent rescaling of F and w: 

F(t, ©) = aF 0), a=r "P, par? 5) 
where D = 2.8. This means that the bath is self-similar 
and can be fully described by a function F(w,). Be- 
low we demonstrate that Eq. (3) is an attractor solution: 
kinetic evolution generically approaches it at large t. 

It is worth noting that Bose star formation can be 
perceived as a second-order critical phenomenon. First, 
its order parameter Mps(t) grows from zero at t > tgr. 
Second, the bath spectrum has thermal small-w tail 
Faw '/? at t ~ tgr, see Fig. 1(b). In Supplemental 
Material A (SM-A) we show that this entails power-law 
field correlators at large distances. The thermal parts 
remain in the self-similar spectra at t > tgr, cf. Fig. 1(c). 

8. Self-similar attractor. Let us ignore the effect of 
Bose star gravitational field on the bath. Then evolution 
of F(t, w) at w > 0 is governed by a homogeneous and 
isotropic kinetic equation [5, 13] 


0,F =StF, (4) 


where St F is the Landau scattering integral — functional 
of F(©) at given 7, see its explicit form in [5] and SM-B. 

Dramatically, the ansatz (3) passes through Eq. (4) at 
any D leaving a one-dimensional equation for the profile, 


(2/D — 1) (w.0y, Fs) — F,/D = St F, , (5) 


This is guaranteed by the scaling St F = a38 St F, re- 
flecting long-range nature of gravitational scattering, 


see [5] and SM-B. The scaling is generic: one can find 
it even using the estimate St F ~ F/tgr. 

On the other hand, Eq. (3) is not a solution if the bath 
is isolated. Indeed, self-similarity gives time-dependent 
mass M, x T*™ and energy Ep x T*£ with 

km =1-3/D, kp =2-5/D, 3kg—5km =1. (6) 
This contradicts to the conservation laws. 

But the ongoing condensation radically changes the 
boundary conditions for the bath. Indeed, the bath 
bosons may scatter, loose energy, and append either to 
the Bose star or to one of its bound states at w < 0. Be- 
sides, with time the star gravitational well grows deeper 
and adiabatically drags low-energy particles to w < 0. 
Both mechanisms absorb bosons with w œ% 0, since grav- 
itational scattering is more effective at low transfers 
Aw & wo [13]. This heats the remaining ensemble due to 
energy conservation. As a result, the bath has decreas- 
ing Mp and growing Ep, i.e. 5/2 < D < 3 in Eq. (6). 

To account for condensation in Eq. (5), we impose a 
condition of finite particle flux at w, + 0 and add an en- 
ergy source StF, > StF; + Js(ws) to the right-hand 
side. This gives a family of solutions at D > 5/2; see 
SM-C for details. The solution F;(w;) with D = 2.8 and 
properly selected J,(w,) is shown in Figs. l(a), (c) by 
chain points. Having almost constant condensation flux 
at low ws, it nevertheless considerably differs from the 
power-law Kolmogorov cascades [14]. 

It is crucial that the self-similar solutions (3) are at- 
tractors of kinetic evolution. This property is apparent 
in Fig. 1(c), but we confirm it explicitly in SM-D by 
solving the full kinetic equation (4) with time-dependent 
source J(r, &). Even if J is essentially non self-similar, 
the solution F(T, @) approaches Eq. (3) with some D. 

4. Growth of the Bose star. In our problem, self- 
similarity of the bath is broken by the Bose star which 
injects energy at its own, non scale-invariant rate J. On 
the other hand, the self-similar solutions are attractors. 
This implies an “adiabatic” regime which was never stud- 
ied before: the bath remains almost self-similar at all 
times, but its parameters slowly drift with time. 

In the first — crude — approximation we can account 
for time dependence of D = D(r). Define km(T) = 
dln M,/dlnr and kg(r) = dlnE,/dlnr. We assume 
that they satisfy the self-similar law (6), 3kg — 5km 7 1, 
if they change slowly. Then the conservation laws M, = 
M — M, and Ey, = E — E, give dint = 3dln(E — E,) — 
5dln(M — M,.) or, integrating, 


(1— E/E — M/M) = (T= T;)/Ta, (7) 


where 7, is an integration constant, M, and E, are the 
parameters of the “dressed” star, and we recalled the 
time translations T > T — 7; [15]. 

To extract the Bose star mass evolution from Eq. (7), 
we estimate the contributions of the excited discrete lev- 


els at w < 0. Theory suggests that large-mass conden- 
sate cannot be accumulated on those levels: it would be- 
come unstable once gravitationally self-bound [16, 17]. 
Then Me < M,;. This is confirmed by our simula- 
tions: Me(t) = Ms — Mos is small and almost constant 
in Fig. 2(a) at r = 2. Moreover, the excited levels with 
w < 0 carry negligibly small energy in simulations as 
compared to the Bose star itself: E, ~ Eps = —y Mẹ, 
where y ~ 0.0542 mG? [17]. Indeed, in Fig. 1(a) only 
the bound states with w œ~ 0 are occupied. Taking 
E, ~ Ers [18] and constant x. = M./M, we obtain the 
growth law for £as(T) = My,/M: 


(1+ a},/e)8(1— ve — Eos)” ~ (T — Ti)/Ta. (8) 


Here e? = E /yM? is a combination of the total mass and 
energy proportional to the invariant E from Refs. [19, 20], 
while M, (T) = (£es + £e) M. 

Note that Ti, T4, and ze in Eq. (8) are empiric fitting 
parameters. However, 7 œ% (1 —7;)(1 — xe)” is fixed by 
the initial condition M,, = 0 at 7r = 1, while ze is small 
and can be ignored, if unknown. This leaves only 7; to 
fit; in fact, 7; ~ —0.1 agrees with all simulations in Fig. 2. 

In Fig. 2(a) we show that the theory (8) (dashed 
lines) reproduces the simulation results for My,(7) and 
M,(r) (solid). A significant statistical test is shown in 
Fig. 2(b) where we display Mp.(7) for 11 simulations with 
e = 0.074 and 22 simulations with e ~ 0.186 (solid data 
vs. dashed theory). These runs have essentially differ- 
ent parameters and kinetic times 10° < wotgr < 3-104. 
Nevertheless, their graphs in Fig. 2(b) merge into two dis- 
tinct curves at two values of e, which agree with Eq. (8). 
Another strong test is performed in SM-E by consider- 
ing self-interacting bosons. In this case our theory still 
describes numerical data, although Eq. (8) gets modified 
by the Bose star self-interaction energy. 

For gravitationally self-bound bath, poL ~ 5/e. This 
means that kinetic approach is valid at €e < 1 [5]. 
Present-day simulations [3, 10] are restricted to € = 0.05, 
cf. Fig. 2. At these values, the Bose star growth is “adi- 
abatic” from the start: dkm,m/dln(T — 7;) < 0.03/e < 1. 
At smaller €, adiabaticity is met at later stages. 

5. Core-halo relation and beyond. At t ~ tgr the en- 
ergy of the baby Bose star is negligible, since Eps x Mj. 
Then Eq. (7) gives M,(t) ~ M (rT — 1)/57, — linear 
growth law for the “dressed” object. 

During longer initial stage, the star is still small, 
Ess < E, and Eq. (8) linearizes to 3x3,/e? + 525s © 
(T — 1)/t.«. Hence, at x, Z € the evolution slows down 
to Mps x t!/3. This transition happens at My, = «eM 
when the virial velocities of the Bose star and the bath 
equalize, |Fy,/Mys| = E/M, i.e. precisely at the “core- 
halo” point of Refs. [21, 22]. The time to the slowdown 
is short at small e: t — ty, ~ Yetg,. This explains, why 
the stars with Mes ~ eM form in cosmological simula- 
tions [3, 21] and seemingly do not grow any further. In 
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FIG. 2. (a) Bose star mass Mp,(t) and the mass of the 
“dressed” object M,(t) in the box simulation of Fig. 1 with 
e ~ 0.074. (b) Evolutions of M,,(t) in 11 + 22 box simu- 
lations with essentially different tg, at two values of e (two 
upper graphs). Numerical results are shown by solid lines, 
while dashed is the theory (8) with ze œ~ 0.043 and 0.021 
at e ~% 0.074 and 0.186, respectively. Circles average over 
simulations with given e. The lower graph shows 2 miniclus- 
ter simulations at e = 0.066 (solid lines) fitted by the theory 
with ze © 0.026 and r; = —0.1 (dashed). For visualization 
purposes, we rescaled the minicluster lines by MM; > 0.7 Mbs. 


truth, the growth continues — hence scatter [19, 23, 24] 
in the simulation results for Mps. 

The next slowdown in Eq. (8) occurs at Ess, = E and 
My, = @/°M. Such heavy objects were observed in 
some cosmological simulations [20, 23]. After this point, 
Mps x t!/. Together, our laws Mp, œ tt/3 and t!/9 agree 
with numerical data in [8] and [6, 7]. 

6. Bose star growth in a halo. Now, consider a denser 
gas which quickly forms a gravitationally bound halo/ 
minicluster under Jeans instability [5, 7, 19]. Using the 
distribution F at w < 0, we find the minicluster mass M, 
its virial energy Eme < 0, and mean particle energy wo = 
—mEmc/M. We also compute its central density p = p(0) 
and potential U(0). This gives tgr, the energy E = Emc— 
U(0)M> 0 counted from the lowest level inside the halo, 
and e? = E/yM?; see details in SM-F. 

With time, the minicluster gives birth to a Bose star. 
The growing mass of the latter is shown in Fig. 2(b) for 
two simulations with e ~ 0.066 (thin solid lines). No- 
tably, M»s(t) is still described by the self-similar the- 
ory with 7; = —0.1 (dashed line), where £e is extracted 
from F. This coincidence strongly supports our theory, 
as it occurs despite the fact that Eq. (8) ignores inhomo- 


geneity of the minicluster. 


7. Discussion. In this Letter we demonstrated that 
kinetics of Bose-Einstein condensation is self-similar if it 
is governed by gravitational (long-range) scattering. This 
solves a long-standing problem [13, 14] with absence of 
Kolmogorov power-law cascades in such systems. The 
Bose star growth law (7), (8) was derived using the new 
working assumption on the “adiabaticity” of scaling ex- 
ponents. This framework may be useful in other contexts. 


To date, simulations of light dark matter structure 
formation [3, 25, 26] cannot provide global distribution 
of Bose stars which are just too small. For that, one 
needs a theoretical input from our Eq. (8) and Refs. [10], 
cf. [27]. Consider, e.g., growing Bose (axion) stars in- 
side QCD axion miniclusters [10]. The latter originate 
from the axion overdensities ® = dp/p|lRp < 10 at the 
radiation-dominated epoch. Equation (8) tells us that 
the star eats the fraction £ẹs ~ 0.1 of the host mini- 
cluster in time t ~ tg,x?,/e®. This time is shorter than 
the the age of the Universe if ® > 6 (10 apsma)?/8@ M34, 
where we used the estimates of [5, 28] and normalized 
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Mis — M/(10-'8 Mo) and mM = m/(10~4 eV) to the 
centers of the discussed minicluster and QCD axion mass 
windows [10, 29]. It is thus realistic to expect that large 
parts of the densest miniclusters are nowadays engulfed 
by their axion stars. Note that the latter may lead to 
spectacular observational effects, see e.g. [30]. 

The other popular model describes growth of gigantic 
Bose stars inside “fuzzy” dark matter galaxy halos [3]. 
Such stars do not reach the “core-halo” point if the re- 
quired time At ~ Qet,, exceeds the age of the Universe. 
This happens if m > 6- 10-24 eV - val? Mg”, where 
we normalized the virial velocity v30 = v/(30 km/s) and 
mass Mg = M/(108 Mo) to the smallest dwarf galax- 
ies. We see that the “fuzzy” Bose stars should be under- 
grown in all galaxies, if the current experimental bound 
m = 2.107% eV [31] on the particle mass is satisfied. 

This Letter is dedicated to the memory of Valery 
Rubakov and Vladimir Zakharov. We thank J. Chan, 
J. Niemeyer, X. Redondo, and S. Sibiryakov for discus- 
sions. The work was supported by the grant RSF 22-12- 
00215 and, in its numerical part, by the “BASIS” foun- 
dation. 


Supplemental material on the article: 
Self-similar growth of Bose stars 


A. DISTRIBUTION FUNCTION 


In the weakly coupled gas, the field Y(t, £) evolves 
almost freely in the mean gravitational field which, in 
turn, changes slowly due to rare scatterings. This means 
that at timescales At < tgr we can write U ~ (U(a)) and 


b(t, 2) ~ X. fn Pala) e. (S1) 


Here {Yn, wn} is the instantaneous eigenspectrum of (U) 
and |fn|? are the occupation numbers of levels wp. Sub- 
stituting (S1) into Eq. (2), we get, 


Fx m$ [fal 6(w — wn), 


where 6(x) = e~(®4+/2)’ At/y4r is the smoothed ô- 
function. This confirms that Eq. (2) defines the distri- 
bution function F ~ dM/dw, indeed. Resolution Aw ~ 
(AtŁ)T} of the latter is of order t7 < Aw « wo. 

In Sec. 2 of the main text we mention that the dis- 
tribution of particles in the box acquires thermal low-w 
tail F xw!/? at t © tgr. Let us show that this entails 
power-law correlator of the field at large distances. Con- 
sider the bath of unbound particles in the periodic box: 


(U) =0, Yn = L73 PaE, wy = p2./2m, Pa = 2nn/L, 
and n € Z3. We assume virialization, i.e. statistical in- 
dependence of different modes within the gas. Then the 
correlator of mode amplitudes equals, 

(fn fj") = Fp) Cae = ne Dd F(wn)/(m? L? pp) , 
where the mean phase-space density f(pn) = (| fanl?) is 


expressed via F ~ (F} which is already time-averaged in 
the definition (2). Equation (S1) gives the field correlator 


(Y(t, wb" (t, y)) = 


27? Í dp F(wp) e'P(@—-y) 
m?L3 J (2r) p ` 


Now, we substitute thermal low-w asymptotic 
F — Fow-'/? at t & tgr, where Fo is proportional 
to the effective temperature. At large |æ — y| this 
corresponds to a power law, 


TFo 


=] 
Vamanga O 2) 


(W(tgr, @)V" (tors Y)) © 


Generically, such power-law behavior is a benchmark of 
second-order critical phenomena. This strongly suggests 
that Bose star formation is a sister process. 


B. LANDAU SCATTERING INTEGRAL 


Let us review Landau kinetic equation for the homoge- 
neous and isotropic gas of gravitating waves [5]. In terms 
of a dimensionless energy distribution F'(r, &), it has the 
form (4), where 


St F = —8;5(r, ©) (S3) 


is the scattering integral related to the Landau flux Š; 
hereafter we mark all dimensionless quantities with 
tildes. The flux  — a cubic functional of F at a 
given T — describes interaction—induced drift of parti- 
cles in the phase space: 


MY (2 2d. BR aoe 
ee fa BA Jnr. (S4) 


3 
Here b is the numerical coefficient from tgr, whereas 
i ee min?/? (©, 0") 
A(@) = f dis’ F?(a') agg (S5) 
Bio) = | da’ F(a’), (S6) 


0 


see Ref. [5] for derivation and details. 

For us, the most important property of the Landau 
scattering integral is its behavior under the scaling (3). 
Substituting the latter into Eqs. (S4), (S5) and chang- 
ing integration variable to &, = Bd’, we find A(@) = 
a”A,(80)/B and B(®) = aB,(6m)/B, where A, and Bs 
denote integrals with F > F,. We get $(@) = a3S,(8a) 
and 


St F(@) = a3 b St F, (B). (S7) 


This last scaling law is used in Sec. 3 of the main text. 

Note that the scaling properties of the scattering in- 
tegral can be understood in a simpler and more gen- 
eral way. To this end we partially restore dimensionful 
units F = MF /2wo, w = 2wo0, and rewrite kinetic equa- 
tion (4) as &F = M St F/(2wotgr). Instead of rescal- 
ing F and & via Eq. (3), we can now change units: 
wo + wo/B and M > aM/8. This gives tgr > tgr / (a°) 
and the same transformation law of the right-hand side 
as in Eq. (S7). 


C. SELF-SIMILAR PROFILES 


In the main text, we introduced two modifications 
of the profile equation (5) to account for condensation. 
First, we impose absorbing boundary condition at w 7% 0: 
enforce F; =0 at ws < wrr and then send the regula- 
tor wrr to zero. We will see that this corresponds to a 
finite and negative particle flux at small w,. Second, we 


FIG. S1. 
Ss(wir) = —1 and (b) D = 2.8, Js(ws) = Jo sech? (ws — w1), 
Jo ~ 0.052, and w; = 1.2. Numbers near the graphs give the 
values of wrr. 


Self-similar profiles with (a) D = 5/2, Js = 0, 


mimic energy income from the condensing particles by 
adding the source J, to the right-hand side of the equa- 
tion, 


(2/D — 1) wsôwu, Fs — F;/D = —ôu, Ss + Js(ws). (S8) 


Here the scattering integral is expressed via Landau 
flux Ss(ws) and the subscript s means that the flux and 
its sub-integrals As (ws), Bs(ws) in Eqs. (54) — (S6) are 
calculated using F, and ws instead of F and ð. 

We turn Eq. (S8) into a set of first-order differential 
equations. First, the definitions (54) — (S6) of the scat- 
tering integrals imply that 0,,,A, = —As/2ws + Cs and 
3v, Bs = F,, where 0,,,C, = —F?/2w,. Second, Eqs. (S3) 
and (S8) can be viewed as expressions for 0,,, fs; and 
u, Ss, respectively. This totals to five equations for the 
unknowns F,,, As, Bs, Cs, and Sz. 

The absorbing boundary conditions imply F; = Bs = 
3A, — 2wiRCs = 0 at ws = wr. They leave two Cauchy 
data C,(wrr) and S,(wrr) which serve as shooting param- 
eters. We tune them to ensure regularity: Fs, Cs > 0 
as Ws > +00. 

At Js = 0 and D = 5/2, the profile equation has a 
scaling symmetry F, > aoFs(ws/ad) with arbitrary ao. 
This is the case when both conditions at ws — +oo can be 
satisfied by choosing C(wrr), while the flux Ss (wr) 4 0 
remains unfixed. If the source is nonzero and D > 5/2, 
the symmetry is absent, and we obtain one solution per 
every D and J;(ws). 

In Fig. Sl(a) we show the solutions F,(ws) with 
D = 5/2 and S;(wrr) = —1, while Fig. S1(b) visualizes 
the case Js 4 0 and D = 2.8. It is clear that the self- 
similar profiles have definite limits wyg — 0. Indeed, 
Eq. (S8) suggests an infrared asymptotic [32] 

F, = Fos w N? + Fis Ws + Oa) ... as Ws > 0 (S9) 
in the unregularized case wy; = 0, where Fos, Fis are 
constants. Imposing this behavior, we obtain the wrr = 0 


FIG. 82. 
Schrédinger-Poisson simulations (pale solid lines). 


Evolutions of the ratio E? /M; in all our 
These 
runs have essentially different tgr, see Sec. 4 of the main 


text. Chain points show the law (7) with r; = —0.1 and 
Ta =1.1-(1—2-)°. 


graphs in Fig. S1 (chain points). Note that Eq. (S9) in- 
cludes a thermal tail at ws — 0, which is indeed observed 
at wir K ws < w in the full numerical simulations, see 
Fig. 1(c) from the main text. 

The profile with wm = 107? from Fig. S1(b) is re- 
peated in Figs. l(a) and (c) of the main text. It has 
Js(ws) = Josech? (ws — w1), Jo ~ 0.052, and wi = 1.2. 
These parameters are selected to fit the simulation data. 

Another good remark is that the self-similar profiles 
fall off as fast as F, x w4 eTo at ws > +oo, where 
q = (4—D)/(2D—4) and ¢ = (2D—4)/ jim (5DA,/w.). 
The cutoff appears because gravitational scattering is in- 
effective at high w. Above the cutoff, the particles cannot 
participate in self-similar dynamics. 

To summarize, the profile equation (S8) has two fam- 
ilies of nontrivial solutions: one solution per every Js at 
D > 5/2 and a branch of D = 5/2 solutions with arbi- 
trary condensation flux S,(wrR) and zero source. 


D. ATTRACTING TO SELF-SIMILAR 
SOLUTIONS 


Let us demonstrate that the self-similar solutions (3) 
are attractors of kinetic evolution. 

To warm up, we explicitly test Eq. (7) using full 
Schrédinger-Poisson simulations. Recall that this law 
approximately describes self-similar bath with slowly- 
varying D = D(r). We compute the bath masses M;(t) 
and energies E»(t) for all our solutions from Figs. 1 and 2 
using the distribution functions F(t, w) at w > 0. This 
gives M, = M — M, and B = E — E,. Then in Fig. S2 
we plot the left-hand side of Eq. (7) versus the right- 
hand side, i.e. basically £3/M? as functions of 7 (pale 
solid lines). The resulting curves are close to each other 
and are well described by Eq. (7) (chain points) despite 


essentially different parameters of the solutions. This 
confirms self-similar character of the bath evolution and, 
hence, our theory for Bose star growth. 

Next, we consider time-dependent kinetic equa- 
tion (4), (S3). To account for condensation onto the Bose 
star, we introduce an absorbing sink W at w ~ 0 and an 
energy source J, 


0,F = -ző —W(a)F + J(r, ©). (S10) 


In practice, we use the sink profile W = Woe @/ wir)” 
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FIG. $3. Numerical solutions of the modified kinetic equa- 
tion (S10). Figure (a) is plotted for J = 0, while (c) con- 
siders self-similar source J = a(r)Js(B(7))/(7 — Ti) with 
Js(ws) = Josech?(ws — wi), Jo ~ 0.1, wi = 1.2, D = 2.8, 
and r; = —0.1. (b), (d) Transformations (3), (S11) of the 
spectra (a) and (c) with parameters (b) D = 5/2, ti ~ —1.9 
and (d) D = 2.8, tri = —0.1 (lines). Chain points show self- 
similar profiles F,(ws) with regulators (b) wrr = 2- 107% and 
(d) wr = 4.5: 1074. (e), (f) Evolutions of D = D(r) and 
E? / MĚ in the case of essentially non self-similar source (lines). 


concentrated at © S wip. It effectively destroys low- 
energy particles at Wọ = 200/wip and wip = 4- 1074. 
We start simulations from the Gaussian-distributed (viri- 
alized) initial state F (0, ©) x &"? e7? and play with 
different J(r, &). 

Figure S3(a) shows the numerical solution F (r, ©) at 
J=0. This is the case when the energy of the bath 
is (almost) conserved and the mass is not: recall that 
the sink swallows particles with w ~ 0. The self-similar 
profile with such properties has D ~ 5/2, see Eq. (6). 
In Fig. S3(b) (dash-dotted and solid lines) we perform 
self-similar rescaling of the spectra (a) with D = 5/2 and 


“UP, Bln) = (r-n), (S11) 


alr) =(r- 7) 
where the time-translation parameter 7; is restored as 
compared to Eq. (3). We see that for properly adjusted 
Ti ~ —1.9 all of the rescaled graphs except for the one 
with 7 = 0 merge into a single curve coinciding with D = 
5/2 self-similar profile F’,(w,) (chain points). It is worth 
noting that the absorbing sink is implemented differently 
in our calculations of F, and F — hence the difference 
in their infrared regulators wig and wip. We see that 
although the starting distribution does not resemble the 
self-similar profile at all, the evolved spectra approach 
aF,;(BÕ©) at T ~ 3 and remain close to it at later times. 
This proves that the self-similar solution with D = 5/2 
is an attractor at J = 0. 

Now, add the energy source with self-similar time de- 
pendence: J = a(r)Js(8(r)@)/(r —7;), where Js(ws) has 
the same form as in Fig. 1; D = 2.8 and 7; = —0.1. 
The respective solution F(T, @) of the kinetic equation 
is visualized in Fig. 53(c). At late times, it exhibits the 
self-similar behavior with D = 2.8. Indeed, the rescaled 
spectra in Fig. S3(d) (lines) coincide at + > 3 with the 
D = 2.8 self-similar profile (chain points). Again, we see 
that the self-similar solutions are attractors. 

Note that the scaling weight D of the solution does not 
always correspond to the time dependence of the external 
source. If the amplitude of J is too large or too small, the 
function F (r, ®©) first attracts to the self-similar profile 
with different D. Later, the weight starts to evolve slowly 
until the source-prescribed value is reached. In such a 
case, the dynamics remains approximately self-similar at 
all times but D slowly drifts with 7. 

The latter situation is illustrated in Figs. S3(e), (£), 
where we consider the source J = Jo 0(r) sech? (© — w1) 
switching on at T ~ 25 as U(r) = [1 + el -7)/A7]-1(7 — 
7) ?/3, where Jo ~ 0.017, wı = 1.2, T; ~ —0.6, T! ~ 24, 
and Ar © 2.3. This time dependence of J(r, ©) explic- 
itly breaks the scaling symmetry of Eq. (S10) making the 
parameter D = (83R—5)/(R—2) in Fig. S3(e) jump from 
D = 5/2 in the beginning of the process to almost 3 in the 
end; to plot the figure, we extracted R = dln Es /dln Mp 
from the numerical evolution. Nonetheless, the combi- 
nation E?/M;? [solid line in Fig. S3(f)] becomes almost 
linear in the late-time region where D = D(r) starts to 
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FIG. 54. Mass evolution of Bose stars in the model with 
nonzero self-coupling at e ~ 0.186 and (a) A = 15, (b) A= —6. 
Thin color lines display results of 9+8 simulations, whereas 
thick dashed is the theory with (a) r; ~ —0.51, x. ~ 0.023; 
(b) 7 ~ 0.14, ze & 0.026. For reference, we repeat the the- 
oretical curve with À = 0 and e ~ 0.186 from Fig. 2(b) (thin 
dash-dotted line). The inset of Fig. (a) shows the theoretical 
curves with À = 0 and À = 15 at larger timescales. 


evolve slowly, again — see the linear fits (chain points). 
This confirms that the solution attracts to self-similarity 
even after the strong kick at T ~% 25. 

It is worth noting, however, that the tilts of the linear 
graphs in Fig. $3(f) are different at early and late times. 
This implies that Eq. (7) holds, but the parameter Tx 
insubstantially changes with time. The latter change is 
ignored in the main text but should be taken into account 
in the refined approaches. 

To summarize, we numerically proved that self-similar 
solutions (3) are attractors of kinetic evolution with a 
sink at w ~ 0 and an energy source. 


E. SELF-INTERACTING BOSONS 


A highly nontrivial test [33] of our theoretical frame- 
work can be performed by studying growth of Bose stars 
in the bath of self-interacting bosons. We describe self- 
interactions by adding the term A|? |Y /(8m?) with cou- 
pling constant A to the right-hand side of the upper 
Eq. (1). This upgrades the full set to Gross—Pitaevskii-— 
Poisson system. 

In the presence of gravity, a comparative effect of self- 
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FIG. S5. (a) The times wotgr of Bose stars formation in mini- 
clusters versus the relevant factor in the theoretical expression 
for this quantity; A = In(poR) is the Coulomb logarithm and 
R is the minicluster size. Filled and empty points correspond 
to simulations of Ref. [5] starting from the 6—distributed gas 
and our solutions with Gaussian initial conditions, respec- 
tively. Thin solid line is the theory with b = 0.7. (b) Time- 
dependent potential U(0) in the minicluster center (dimen- 
sionless units). Thin solid graphs are extracted from two long 
simulations, while the dashed line is the time-averaged value. 


interactions is characterized by a dimensionless combina- 
tion [17] A = 2Awo/(m3G), where wo is the typical par- 
ticle energy. Our self-similar solution (7) is applicable if 
gravity dominates, i.e. at [5, 34] 


Tor or ? 


~ ~ 1. 
TA Ogr 10247? In(poL) = 


(S12) 


Here Tgr and 7) are the gravitational and self-interaction 
relaxation times, while og, and a) are the respective 
transport cross sections. Below we keep [35] Tg,/T, < 
107? in all simulations. 

Despite Eq. (S12), self-interactions can modify the 
growth law of Bose stars, and their effect can be consid- 
erable [7]. Indeed, although the terms proportional to A 
are small inside the light stars, they grow with mass and 
start to dominate [36] at My, > M, = |AG|~1/?. If the 
self-interactions are repulsive, Aà > 0, this increases the 
energy of heavy stars to Eps x —M?,. The attractive case 
is more dramatic: negative self-pressure makes the Bose 
stars collapse [7, 36-38] as Bosenovas at My, Z, 10.2 My. 
Generically, we write 


Ess = —yM3, €(Mys/My), (S13) 


where the function € accounts for self-interaction energy. 
It is smaller (larger) than 1 at A > 0 (A < 0). Besides, 
E x1 at Mes < My when the self-interactions are negli- 
gible. In practice, we compute € numerically by solving 
the Gross—Pitaevskii—Poisson system for every Mp,/M). 

Using Eq. (7), we obtain the growth law of self- 
interacting stars [cf. Eq. (8)], 


(1+ REP (1 — te — Los) & (T — Ti)/Ta, (S14) 


where E = E(a,,M/M)). At the qualitative level, 
Eq. (S14) agrees with the phenomenon suggested in 
Ref. [7]: for fixed 7, and 7; the Bose stars grow faster 
for positive A (£ < 1) and slower for negative \ (E < 1). 
This feature is illustrated in Figs. S4(a), (b) that com- 
pare two theoretical curves £ps(T), Eq. (S14), at A = 15 
and A = —6 (thick dashed lines) with the one at A = 0 
(thin dash-dotted). 

We performed an explicit numerical test of Eq. (S14). 
Starting from the Gaussian-distributed gas in the box, 
we performed 9 simulations at À = 15 and 8 simu- 
lations at à = —6. We used M = 20po/m?G and 
L = (35 + 40)/po. Mass evolutions of the respective Bose 
stars are shown in Figs. S4a, b by thin color lines. They 
are well described by Eq. (S14) (thick dashed). How- 
ever, the values of the fitting parameter [39] 7; are differ- 
ent with respect to non self-interacting case: we obtain 
T; œ% —0.51 at \=15 and 7; ~ 0.14 at À = —6. 

It is worth noting that the dependence of 7; on À affects 
the growth law of Bose stars. At moderately small 7, it 
may even compensate the (de)acceleration effect of self- 
interaction energy, see the inset in Fig. S4(a). However, 
at large timescales the self-energy wins and makes the 
growth go faster at A > 0 and slower at A < 0 — see the 
inset, again. 


F. SIMULATIONS IN MINICLUSTERS 


Although the application of our theory (8) is straight- 
forward at the qualitative level, things become more 
tricky once precise agreement with simulations is re- 
quired. To this end, we accurately determine the mini- 
cluster parameters. 

We form gravitationally bound miniclusters by trig- 
gering strong Jeans instability in the dense virialized 
gas [5, 7]. In particular, our two long simulations start 
from very large mass Myo = 112.5/wo in the box L = 
52.5/po. At these values, the miniclusters engulf more 
than 55% of matter, and the remaining diffuse particles 
do not affect much the growth of objects within them. 

We define the minicluster center as the center-of-mass 
of matter distribution within the box; we call it æ = 0 for 
simplicity. The density p = (0) in the minicluster center 
is then obtained as the value of p = m|w(t, x)|? averaged 
over the Gaussian spatial window. The remaining pa- 
rameters are extracted from the distribution function (2) 
or, specifically, from its part at w < 0 that describes a 
self-bound minicluster. Namely, the mass M and energy 
Eme < 0 of the minicluster are obtained by integrating 
F and wF/m over this region. Then the virial parti- 
cle energy equals wọ = —mEmc/M, the virial radius is 
R = (3wo/2mmGp)'/?, while A = In(poR) is the Coulomb 
logarithm. 

Once the minicluster parameters are specified, we find 
the numerical factor [40] b in the expression for the relax- 


ation time tgr. To this end we perform many short-time 
simulations at different values of parameters and wait 
until Bose stars appear in their miniclusters. The mo- 
ments t,, when they form (empty points in Fig. $5(a)) 
are well described by the theory with b = 0.7 (line) — 
the same value as in Ref. [5]. The coincidence of b’s 
is remarkable because minicluster simulations of Ref. [5] 
(filled points) start from the d-distributed gas in the box, 
|p|? x 6(|p| — po), while our simulations use Gaussian 
gas with |p|? œx e~?’/0. This suggests that formation 
of miniclusters strongly intermixes the gas forcing it to 
“forget” the initial condition. 


An important part of our procedure is a computation of 
the gravitational potential U(0) in the minicluster center. 
In the notations of Eqs. (7), (8), this parameter enters 
the total energy E = Eme — U(0)M which is positive 
and counted from the lowest level inside the minicluster. 
Notably, the value of U(0) visibly drifts with time, since 
the minicluster gets eaten by the Bose star and becomes 
lighter; see Fig. S5(b). We calculate the potential using 
the Bose star itself as a sensor. On the one hand, its 
mass Mps and binding energy wbs = —3myM?, can be 
extracted from the profile |~»,(a)|?. On the other, the 
“Bose star” peak in the energy distribution is located 
at w = wbs + mU(0). Subtracting these quantities, we 
obtain the solid lines in Fig. S5(b) corresponding to two 
long simulations. We use the time-averaged value of U(0) 
(dashed horizontal line) in the theoretical expressions for 
E and ê = E/yM?°. 

Finally, we determine the fraction £e = Me/M of par- 
ticles on the discrete levels of the Bose star potential in 
the same way as before: by integrating F over the region 
wbs + mU(0) < w < mU(0). Once this is done, the theo- 
retical predictions (S14) match the Bose star mass curves 
extracted from the simulations (lower graph in Fig. 2(b)). 
The respective best-fit value 7; ~% —0.1 matches that in 
the box simulations. 
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